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ABSTRACT 

We investigate the spatial distribution of charged particles accelerated by non- 
relativistic oblique fast collisionless shocks using three-dimensional test-particle sim¬ 
ulations. We find that the density of low-energy particles exhibit a localised enhance¬ 
ment at the shock, resembling the “spike” measured at interplanetary shocks. In 
contrast to previous results based on numerical solutions to the focused transport 
equation, we find a shock spike for any magnetic obliquity, from quasi-perpendicular 
to parallel. We compare the pitch-angle distribution with respect to the local mag¬ 
netic field and the momentum distribution far downstream and very near the shock 
within the spike; our findings are compatible with predictions from the scatter-free 
shock drift acceleration (SDA) limit in these regions. The enhancement of low-energy 
particles measured by Voyager 1 at solar termination shock is comparable with our 
profiles. Our simulations allow for predictions of supra-thermal protons at interplane¬ 
tary shocks within ten solar radii to be tested by Solar Probe Mission. They also have 
implications for the interpretation of ions accelerated at supernova remnant shocks. 

Key words: Physical Data and Processes: turbulence; ISM: cosmic rays, magnetic 
fields. 


1 INTRODUCTION 

Collisionless shocks are considered major sources of ener¬ 
getic particles in interplanetary and interstellar medium. 
The observational association of energetic particles with 
shocks ranges from (keV-MeV) for solar energetic particles 
originating in solar events (Lario et al. 2003; Giacalone 2012) 
to the galactic cosmic rays whose TeV signature has been 
observed in extended shell-type supernova remnants (Albert 
et al. 2007; Acciari et al. 2011). 

Localised enhancements of energetic particles are often 
seen at the time of passage of interplanetary shocks by en¬ 
ergetic particles detectors. As recently reported by Lario et 
al. (2003), in addition to the classic energetic storm parti¬ 
cle events (a few hours), structures seen at travelling inter¬ 
planetary shocks, e.g., by ACE energetic particle detector, 
also spike events (~ 10 min or shorter) are observed. En¬ 
hancements of energetic proton intensities at shocks were 
observed by Explorer 33-35 (Armstrong et al. 1970) or Vela 
4 (Singer & Montgomery 1971). Axford & Reid (1963) inter¬ 
preted these as evidence of anisotropic reflection of energetic 
particles between Earth’s bow shock and shocks travelling in 
the solar wind. This interpretation was questioned by Arm¬ 
strong et al. (1970) who found in some events energetic- 
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particle features behind the shock similar to the upstream. 
Spike events were modelled by Decker (1983) who considered 
shock drift acceleration (SDA) of a seed particle population 
at a quasi-perpendicular shocks (Obu > 70°, where Obu is 
the angle between the shock normal and the upstream or¬ 
dered magnetic field). However, an extension of the SDA sce¬ 
nario in Decker (1983) was needed to explain the energetic 
storm particle events, or more generic post-shock enhance¬ 
ments, associated with Obu < 70°. Scholer (1985) argued 
that the self-excited waves produced upstream by energetic 
ions via streaming instability might account for the local en¬ 
hancements at shocks with Obu < 70°. Upstream magnetic 
traps for energetic particles formed by multiple encounters of 
the field lines with the shock surface have also been invoked 
as a way to generate spikes through the collapse of the trap 
into the shock as the field advects with the fluid (Erdos & 
Balogh 1994). The conditions on the shock obliquity under 
which particle anisotropy leads to shock spike were analysed 
in Gieseler et al. (1999). The crossing of the solar termina¬ 
tion shock by Voyager-1 (Decker et al. 2005) also seemed to 
exhibit a spike in high-energy protons (3.4— 17.6 MeV) and 
ions (40 — 53 keV) intensities, as was suggested by le Roux 
et al. (2007) and Zuo & Eeng (2013). 

The spatial scale of the increase of particle intensity 
that we refer to as a “spike” is of the order of gyroscale, sug¬ 
gesting that the diffusion approximation does not apply to 


2 F. Fraschetti and J. Giacalone 


the spike formation. Thus, there must exist an anisotropy in 
the pitch-angle distribution. The origin of energetic particles 
intensity spikes across collisionless oblique non-relativistic 
shocks have been considered in several recent studies, using 
two different approaches with diverging outcomes: 1) Monte- 
Carlo simulations (Ellison et al. 1995), i.e., the trajectories of 
a large number of particles are simulated with an ad hoc pre¬ 
scription for the pitch-angle scattering (Ellison et al. 1996) 
with no emerging spike; 2) focused transport equation, or 
ETE, (Isenberg 1997; Kota & Jokipii 2004) which relaxes the 
assumption of near isotropy in the pitch-angle distribution 
and does not assume continuity of the particle distribution 
across shocks (see, e.g. le Roux et al. (2007)). Contradictory 
results about the formation of shock spike were reported also 
for relativistic shocks, as shown by Monte Carlo test-particle 
simulations (Ostrowski 1991; Naito & Takahara 1995). 

Observational consequences of spikes are relevant 
also to astrophysical shocks. Piling-up of energetic elec¬ 
trons at the shock might result in an enhanced filamen¬ 
tary synchrotron-emitting structures at supernova remnant 
shocks, as argued in Gieseler et al. (1999). An enhancement 
of particle density at the shock might also affect filamenta- 
tion structures thought to be seeded by Weibel instability 
at relativistic shocks (Medvedev & Loeb 1999). 

The spike-like discontinuity found in the upstream par¬ 
ticle density was interpreted as due to particles reflected 
after one shock-crossing and conserving the magnetic mo¬ 
ment as they encounter the magnetic barrier at quasi¬ 
perpendicular shocks (Gieseler et al. 1999). However, for 
smaller obliquity shocks this interpretation does not seem to 
apply. Previous approaches do not clarify the dependence of 
the spike intensity on Obu- A shock spike emerges from the 
solution to the ETE regardless of Obu] however, for quasi¬ 
parallel shocks it has been argued (le Roux et al. 2007) that 
the spike is unphysical and it should disappear if the par¬ 
ticle momentum were calculated in the shock rest frame, 
instead of being calculated in the local plasma frame, as in 
the Parker (1965) approximation. 

In this paper we present results from the numerical in¬ 
tegration of test-particle trajectories in kinematically pre¬ 
scribed turbulent fields associated with a collisionless shock. 
The geometry comprises a planar shock with an arbitrary 
obliquity of the upstream average field. Particle scatter¬ 
ing and decorrelation arises from directly integrating the 
particle equation of motion in the given turbulent back¬ 
ground, without ad hoc prescription of pitch-angle change. 
This procedure provides us with spatial profiles of the in¬ 
tensity of energetic particles with an unexpectedly signifi¬ 
cant enhancement at quasi-parallel shock. We compare the 
pitch-angle anisotropy and the momentum distribution at 
the shock with far downstream. Our density profiles are con¬ 
trasted with the DSA prediction in various particle momen¬ 
tum range, from supra-thermal up to energetic particles. As 
in our previous works (Eraschetti & Giacalone 2012), parallel 
and perpendicular diffusion originate from a fully specified 
magnetic turbulence. 

This paper is organized as follows: in Sect. 2 previous 
extensions to non-diffusive regime are outlined; in Sect. 3, 
complemented by the appendix, we describe details of the 
numerical approach. In Sect. 4 we present sample charged 
particle trajectories at the shock, identifying the role of 
magnetic obliquity and turbulence power in the energiza¬ 


tion process. In Sect. 5 we show that spatial density profiles 
exhibit an enhancement at the shock regardless of the shock 
obliquity and for sufficiently low-energy particles. In Sect. 6 
we compare the distribution of the pitch-angle with respect 
to the local fluctuating field ahead and behind the shock 
with the expected near-isotropy assumed by DSA and with 
the model for scatter-free SDA. In Sect. 7 we compare the 
momentum spectrum at the spike with the far downstream 
high-momentum power-law. In Sect. 8 we discuss the im¬ 
plications of our findings for the in-situ measurements of 
interplanetary shocks. In Sect. 9 we discuss the findings and 
outline the conclusions. 


2 NON-DIFFUSIVE APPROACH 

Standard DSA theory associates to shocks a population 
of energetic particles with a pitch-angle distribution nearly 
isotropic in the local plasma rest frame. This implicitly as¬ 
sumes that the speed v of those particles in the local plasma 
frame is much greater than the speed of the shock in the 
upstream fluid Ui. In other words particles have time to 
scatter on the magnetic irregularities and isotropize their 
pitch-angle distribution as the bulk fluid advection occurs 
on much greater time-scale. In the limit that particles ad¬ 
here to magnetic field lines, v is to be compared with the 
shock speed in the upstream fluid projected along the up¬ 
stream average field, i.e. Ui! cosObu^ namely the speed along 
the shock of the intersection point field line/shock. Clearly 
this ought not to be valid if decorrelation from field line is 
allowed. 

The steady-state solution of the Parker transport equa¬ 
tion (Parker 1958) for a one-dimensional planar shock and 
uniform plasma flows on either side predicts that the par¬ 
ticle density rises exponentially from far upstream toward 
the shock, with an e-folding distance depending on the spa¬ 
tial coefficient diffusion n in the direction perpendicular to 
the local shock surface and on the fluid advection speed 
Ui. The intensity is continuous at the shock along with its 
Lagrangian derivative and in the downstream flow it is uni¬ 
formly equal to a constant depending only on the parti¬ 
cle momentum. In the presence of self-generated waves the 
steady state upstream profile rarefies more slowly than ex¬ 
ponentially (l/x) (Bell 1978). 

Deviations from the assumed near-isotropy have been 
already studied in several cases even in the absence of 
shocks. For instance, scattering on small-scale turbulence is 
found to generate an inertial force due to the velocity shear 
of second order in Ui/v (Earl et al. 1988), usually neglected 
with respect to the adiabatic energy change in the trans¬ 
port equation (V • U). If the random motion of scattering 
centres is embodied in the momentum diffusion coefficient, 
the standard second-order \nUi/v Fermi acceleration arises. 
Neglecting orders higher than first in Ui/v requires disre¬ 
garding the anisotropic part, in the local plasma frame, of 
the phase-space particle distribution function, although one 
can solve for the isotropic part of / including second order 
momentum diffusion. 
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3 NUMERICAL SET-UP 

We carry out numerical simulations of supra-thermal ^ test- 
particles (protons) at collisionless fast shocks travelling in a 
medium with an embedded turbulent magnetic field with 
a moderately high upstream Alfven Mach number, i.e., 
Ma ^10 (compression of transverse components of U and 
B depends on Ma)- We consider the limit of the initial par¬ 
ticle speed comparable to the bulk plasma flow speed. In 
our model, particle scattering arises from the interaction 
between the particles and a kinematically prescribed turbu¬ 
lent magnetic field. We examine steady-state particle den¬ 
sity profiles across the shock for various particle energies, 
variance in the turbulent random magnetic field, and Osn- 

The geometry used is depicted in Fig. 1, upper panel. 
The shock surface is at the plane x = 0. Bulk plasma flows 
into the shock from x < 0 to x > 0 along the shock normal 
with speed Ui <C c, measured in the shock frame. Mass 
conservation across the shock implies a decrease of the bulk 
plasma speed from upstream to downstream along XI 0 U 2 — 
Ui/r, where r is the density compression at the shock. The 
profile of the bulk velocity is taken to be discontinuous at 
the gyroscale of the particles considered, much larger than 
the shock thickness, of order of ion skin depth; thus we use 
Ux{x) = Ui for a: < 0, Ux{x) = f /2 for x > 0 (see Appendix). 

The upstream three-dimensional magnetic field is given 
by B(x) = Bq + (5B(x), with an average component Bo hav¬ 
ing orientation Obu with respect to the shock normal and a 
random component ^B = (5B(x,y, z) having a zero mean 
(((5B(x)) = 0); the correlation length of the turbulence, Lc, 
is much greater than the particle gyroradius, r^, at all ener¬ 
gies considered in this paper. 

For the sake of simplicity, we assume that the three- 
dimensional power-spectrum of (5B in the inertial range, 
G{k)^ is isotropic and scale-invariant, or Kolmogorov: 
G{k) oc where k is the wavenumber magnitude, 

jS = 5/3 is the one-dimensional power-law Kolmogorov in¬ 
dex and the additional 2 accounts for the dimensionality of 
the turbulence (see Appendix). 

Non-thermal seed particles are injected at the shock 
with the same momentum magnitude, po = myoUo, and 
pitch-angle isotropic momentum distribution in the local 
plasma rest frame, where the cosinus of the pitch-angle with 
respect to the average field is defined as Ugiobai = P'^o/pBo. 
A sample particle trajectory is shown in Fig. 1, lower panel. 
We solve numerically the Lorentz force equation for the 
charged particles in the prescribed magnetic field (see Ap¬ 
pendix for details). For any Obu, the particle injection speed 
in the local plasma frame vq is comparable to the flow speed 
(uo = 3.5f/i). The analytic solution of the induction equa¬ 
tion is used in the particle equation of motion that is numer¬ 
ically solved with the procedure described in the Appendix. 
As usual, since the plasma is infinitely conductive, i.e., in the 
plasma rest frame both upstream and downstream the elec¬ 
tric field vanish, particles undergo acceleration by the shock 
frame electric field E^(x, y, z, t) = —U(x)/c x B^(x, y, z, t), 
where U(x) = {Ux,0,Uz) and T indicates the plane or¬ 
thogonal to the shock frame flow velocity U(x); in Ap- 


^ Here supra-thermal particles have speed v greater but compa¬ 
rable with the upstream flow speed in the shock frame, Fi, in 
contrast with high-energy particles (v ^ Ui). 



Figure 1. Upper panel - Shock frame configuration of MHD vari¬ 
ables (adapted from Giacalone (2005)). Lower panel - 3D trajec¬ 
tory of a sample particle injected at an oblique shock, propagating 
along x-axis, with Obu = 83° and = 0.3. Spatial axis are in 
units of initial upstream gyroradius r®. Color scale indicate in¬ 
crease of kinetic energy in the local plasma frame. 


pendix B(x,y,z,t) is related to 'B(xo,y, z,to) = Bo(xo) + 
(5B(xo, y, to). We neglect the electric field of the order of 
vaBI c arising from the magnetic field fluctuations modelled 
as Alfven waves propagating along the field at speed ua, 
as measured in the local plasma frame (see Schlickeiser & 
Vainio (1999) for a review on the limits of the magnetostatic 
approximation). This is justified by the fact that for the pa¬ 
rameters used in our test-particle simulations the injection 
particle speed is much greater than Alfven speed on both 
sides of the shock (see Sect. 5 for details). In this work we 
focus on the particle acceleration resulting from the fluid 
compression at the shock. Stochastic acceleration, relevant 
as the particles move away from the shock, is expected to 
have a smaller effect on the spike formation because 1) the 
electric field vaB jc is disordered and 2) the acceleration is 
efficient only in resonant condition, i.e., Xg comparable with 
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0 ^ = 1.0 



Figure 2. Momentum in units of injection momentum po versus 
the ar-coordinate for three particles injected with random veloc¬ 
ity at the shock, located at a: = 0, for various shock obliquities 
6Bn = 1.0). Length is in units of initial upstream particle 

gyroradius; notice the difference in spatial scale for various Osn- 
Upstream medium is to the left. Orbits show multiple shock cross¬ 
ing. 


length-scale of the turbulence, not reached in our simula¬ 
tions (the study of the stochastic acceleration is deferred to 
a separate work). 

In order to obtain a steady-state solution, we use an ap¬ 
proach similar to that in Giacalone (2005) (see Appendix). 
We do not assume a free-escape boundary upstream but, 
instead, follow particles until either a) they reach a pre¬ 
specified high-energy cut-off ot h) escape by advection down¬ 
stream. Energetic particles, once in the upstream medium, 
undergo scattering on the pre-existing magnetic turbulence 
which advects with the upstream plasma. Our assumed mag¬ 
netic power spectrum includes only the pre-existing fluctua¬ 
tions and not those that are self-generated by the particles. 

The sought spatial structure of particle density across 
the shock has a scale of hundreds of supra-thermal particle 
gyroradii. Thus, the test-particle approach can easily cap¬ 
ture the proton scales we are interested in; also, the thickness 
of the shock dissipation layer (of order of ion skin depth) 
is negligible with respect to the gyroradius of the supra- 
thermal particles. 

The relevant parameters are the ratio of the parti¬ 
cle gyroradius to the turbulence correlations length, r^/Tc, 
the normalized upstream turbulence energy density, — 
^ the upstream magnetic obliquity Obu- Therefore 
our treatment applies to energetic particles at astrophysical 
shocks in various contexts: from the interplanetary, to the 
supernova remnant, to the intergalactic medium. 


4 PARTICLE TRAJECTORIES 

In the present section, we show individual particle trajec¬ 
tories to illustrate salient features of the particle energiza¬ 
tion at shocks with various obliquity; in the following sec¬ 
tion we will describe the resulting spatial profile of a mono- 
energetically injected population of supra-thermal particles, 
with random pitch-angle and phase. 

We have numerically integrated the proton trajectories 
for various Obu and . Figure 2 traces the actual trajecto¬ 
ries of sample protons injected at shock for three different 
Obu and <j^ — 1.0. The injection particle speed in the local 
(downstream) plasma frame vq is only 3.5 times the speed of 
the shock in the upstream plasma frame Ui for every Obu- 
Particles are energised preferentially along the shock surface, 
regardless of the value oi 0Bn ^ suggesting that drift accelera¬ 
tion plays a significant role even at shocks with = 0°. If 
the pre-existing upstream magnetic fluctuation is sufficiently 
strong (cr^ = 1.0, Fig. 2), making (5B spatially isotropic, the 
turbulent motional electric field —U/c x (5B has a compo¬ 
nent also along the shock surface and contributes to acceler¬ 
ate particles even at OBn = 0°. If SB is weaker, the particle 
spends more time far from the shock: the acceleration time is 
given by race ~ i^xjUi (Drury 1983), where Kx is the spatial 
diffusion coefficient along the direction normal to the shock. 
From quasi-perpendicular to parallel shock (the latter re¬ 
quiring longer elapsed time than the former as /^y), 

the particle excursion far upstream of the shock increases by 
two orders of magnitude due to the advection of the mag¬ 
netic field lines dragging the particle back to the shock for 
Osn > 45°. 

As for small momenta (p < lOpo) as seen in Fig. 2, 
a parallel shock confines more efficiently particles around 
the shock through an enhanced scattering if cr^ = 1.0 than 
= 0.3. In the next section we confirm that for = 1.0 
the spike ai a OBn > 0° shock is higher (relative to the 
downstream background) and persists at high momentum 
compared to cr^ = 0.3: high cr^ increases the permanence 
time of particles at shock, thus increasing the contribution 
of the shock crossing to the particle acceleration. Conversely 
for high-obliquity shocks: the permanence time within the 
spike increases for smaller as spike is generated by the 
SDA nearly scatter-free process (cfr. Fig. 3), as first mod¬ 
elled in Decker (1983): particle slides along y — z coordinates 
under the turbulent electric field —U/c x B (B = Bo + 5B) 
while remaining close to the shock. 


5 SUPRA-THERMAL PARTICLES DENSITY 
PROFILES 

In this section we present the steady-state spatial density 
profiles N{x) for accelerated protons from our numerical 
simulations. We show N{x) at various momenta as measured 
in the local plasma frame to be compared directly to the so¬ 
lution of the steady state transport equation. Particles are 
injected with volUi = 3.5. A strong shock is used (plasma 
density compression is r = 3.5) with Alfven Mach-number 
Ma = 7. The downstream boundary Xh — 7.5 x 10^r° 
(see Appendix) is chosen to ensure that the highest-energy 
particles have enough space in the downstream medium to 
isotropize in pitch-angle. As a momentum boundary (de- 
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Figure 3. Steady-state spatial density profiles of accelerated par¬ 
ticles for various momentum intervals, and correspondent proton 
energy, at = 85° and cr^ = 0.3. The density is normalised 
to the far downstream value. The distance is in units of initial 
upstream gyroradius at 1 AU. 

scribed earlier), we use pb — 500po- For a strong shock trav¬ 
elling at 1 AU in the solar wind, the ratio uo/Ui = 3.5 
corresponds to a proton with kinetic energy 10 keV at a 
shock propagating with Ui — 400 km s“^ in the upstream 
plasma frame with an embedded background magnetic field 
{Bi) = 5 nT (r° = 2,900 km and Uo = e{Bi)lmc = 0.48 
s“^); the upstream magnetic fluctuation has a correlation 
length Lc = 0.01 AU. Note that the momentum boundary 
Pb is chosen so that highest-energy particle gyroradius in the 
upstream is smaller than the Lc. 

Figures 3, 4 show numerically computed spatial density 
profiles of accelerated particles across the shock for various 
plasma-frame momenta. Unlike the prediction of the stan¬ 
dard diffusion theory, we find a local enhancement (even 
more than two-fold) at the shock for low-energy accelerated 
particles at both quasi-perpendicular shocks (Obu = 85°, 
Fig. 3) and parallel shocks (Obu = 0°, Fig. 4). Note that the 
quasi-perpendicular shock profiles are shown for weak tur¬ 
bulence (cr^ = 0.3) and the parallel shock profiles are shown 
for strong turbulence (cr^ = 1-0), respectively in Fig.s 3, 4; 
those values maximise the amplitude of the shock spike 
(see also Sect. 6). The amplitude of the spike decreases for 
increasing particle energy in both cases, as expected in the 
standard DSA theory for sufficiently fast particles: at higher 
momentum (v ^ Ui) the spike vanishes and the uniform 
downstream density is recovered. The length-scale of the 
intensity decay ahead of the shock depends on the spatial 


Figure 4. Steady-state spatial profiles for various particle mo¬ 
mentum intervals at Ob^ — 9° cr^ = 1.0 (see Fig. 3). 



Figure 5. Density enhancement at the shock in units of far down¬ 
stream value as a function of 6Bn for various cr^ and plasma-frame 
particle momentum. Note the difference in scale of the vertical 
axis. 


diffusion coefficient along x-axis (for the same Ui) thus ex¬ 
tends less for the quasi-perpendicular than the parallel shock 
(typically A^y). However, the presence of a shock spike 

might modify the exponential decay upstream: the compar¬ 
ison of the density profile upstream with the diffusion pre¬ 
diction and local non-diffusive behaviour is deferred to a 
separate work. 
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Our simulations contrast with results from other anal¬ 
yses (see, e.g., le Roux et ah (2007)) which suggest that 
only near-perpendicular shocks should exhibit a shock spike 
and the spike at parallel shock would be an artifact of the 
FTE method. In quasi-perpendicular case our downstream 
tail of supra-thermal particles density arises from the quasi 
scatter-free regime (cr^ = 0.3), as shown in the anisotropy 
analysis in Sect. 6, whereas in the parallel shock it results 
from isotropic scattering (cr^ = 1.0). An additional contrast 
is the qualitative structure of the spike: we find a trapping 
tail in the downstream density with thickness of order 10^r°, 
whereas FTE solutions or Monte-Carlo simulations sharply 
drop from the shock spike to the far downstream limit. The 
thickness of the downstream tail increases from perpendic¬ 
ular to parallel case, as a consequence of k,± k,\\. We con¬ 

clude that the deviation from DSA found here has different 
physical origin from the one discussed in previous works. 
Such a refined shock structure is observable by future inter¬ 
planetary shock missions as discussed in Sect. 8. 

Spike enhancements for various Obu, cr^ and particle 
momenta are summarised in Fig. 5. At high-obliquity shocks 
(Obu ^ 70°), the shock spike is much more enhanced at weak 
than strong turbulence, for every momentum, suggesting 
that the weaker is the fluctuating component of B the higher 
is the shock enhancement, that originates predominantly in 
a scatter-free drift acceleration limit, i.e., for uniform and 
static electric and magnetic field. Such an interpretation is 
corroborated by the pitch-angle anisotropy in the following 
section. The drop in enhancement between Obu = 70° and 
85° can be attributed to a more efficient advection of par¬ 
ticles away from the shock as Obu approaches 90° than at 
Obu = 70°. Conversely, at parallel shock the spike ampli¬ 
tude is increased by strong turbulence. Despite the different 
particle acceleration processes dominating at different Obu, 
a downstream spike emerges at the same scale (~ 10^r°). 
For large Obu, the relative importance of the SDA regime 
decreases as is increased. For small Obu scattering on 
the turbulence on both sides of the shock is the main ac¬ 
celerating agent as maintains particles in the shock vicinity; 
thus the spike is enhanced at large . We therefore argue 
that the FTE solution should find a more prominent spike 
at stronger turbulence (cr^ approaching 1). 

Our steady-state density profiles are not qualitatively 
affected by the shock strength. We performed simulations 
with density compression at the shock smaller than 3.5 (2.5 
and 1.5) and found no qualitative difference in the spike 
enhancement; on the other hand, smaller compression leads 
to softer downstream momentum distribution (see Sect. 7), 
as expected from DSA, needing a larger particle pool. 


6 ANISOTROPY 

In this Section we compare at various locations the plasma- 
frame pitch-angle anisotropy by using a pitch-angle defined 
with respect to the local magnetic field at the particle po¬ 
sition, with no spatial average: fiiocai = P • B/pS. The dis¬ 
tribution f{fiiocai) is depicted in four spatial regions in each 
row of Fig. 6: two shells of thickness 30r° on both sides of the 


Ben = 85°, [2-4] Pq, UPSTREAM 



Figure 7. Synopsis of the steady-state pitch-angle iHocal distri¬ 
bution in the upstream plasma at various shells of 30r® thickness 
at increasing distance from the shock, downward from [—30 : 0]r^ 
as far as [—450 : —420]r®. Two different turbulence strengths 
are compared (cr^ = 0.3 dotted line, cr^ = 1 solid line). Here 
^Bn = 85°, vq/Ui = 3.5 and local plasma frame momentum in 
the range [2 — 4]po- 

shock, far upstream and downstream^. In Fig. 6 we consider 
two distinct momentum intervals ([2 — 4], [10 — 20]po), three 
magnetic obliquities (Obu — 85°,45°,0°) and take — 1.0. 
The far downstream distribution is calculated in the inter¬ 
val between O.lx?, and which is a region much larger 
than the particle mean free paths so that the distribution 
is isotropic. The far upstream distribution is calculated in 
X < — 30r° to seek excess of reflected particles, whose pile-up 
was argued by Gieseler et al. (1999) to originate the spike, 
that they found located solely in the upstream region. Such 
a thickness is broad enough to account for various particle 
populations: low-energy particles drifting along the shock in 
a multiple encounter^ interaction, accelerated particles ad- 
vected downstream, particles coming back to the shock from 
far downstream, particles reflected one or multiple times by 


^ For a proton injected at kinetic energy F = 10 keV in the solar 
wind at 1 AU, the shell has thickness 30r® ~ 6 x 10“"^ AU. 

^ It is customary to define single encounter with a shock the 
period of time during which the particle remains within one gy- 
roradius from the shock. 
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Figure 6. Comparison of the steady-state pitch-angle fJ^iocal distribution in the plasma frame in the far upstream (downstream) flow 
within two shells of thickness 30r® upstream (downstream) of the shock for charged particle propagating in a strong turbulence (cr^ = 1.0) 
at various magnetic obliquity (Obu = 85°,45°,0° in various rows). Here vq/Ui = 3.5. The panels on the left (right) column refer to the 
upstream (downstream) fluid. In every panel the red (blue) curve corresponds to the fluid at the spike (far from the shock). Two ranges of 
momentum are considered ([2 — 4], [10 — 20] in units of po) to show the progressive flattening of the excess at p^iocal = “1 (upstream) and 
of the “pancake” at liiocal = 0 (downstream) for increasing energy. At high momentum, p > 20po5 fhe pancake disappears and f(fJ>iocal) 
collapses to an isotropic distribution, as expected from DSA. The relatively large numerical fluctuations at 6Bn = 85° upstream (top 
left panel) originate from the upstream lower statistics for a quasi-perpendicular shock. 


the shock moving to the upstream and particles advected 
back to the shock from upstream. 

We find that in the downstream shell (Figs. 6, right col¬ 
umn, all panels, red curves) the lower momentum ([2 —4]po) 
pitch-angle distribution has a “pancake” shape in logarith¬ 
mic scale (broad bump peaking at uiocai = 0, i.e., for parti¬ 
cles moving in the direction perpendicular to the local mag¬ 
netic field), whereas in the far downstream fluid f{niocai) 
is closer to isotropic, as predicted by DSA. Likewise, within 
the scatter-free limit (cr^ = 0) at quasi-perpendicular shocks 
a distribution peaked at Ugiobai = 0 emerges (see the fi- 
histograms after a single shock encounter in Figs. 7-17 
of Decker (1988)). Such a peak is broadened by small- 
amplitude perturbations = 0.09 (Decker 1988) approach¬ 
ing the pancake-shape we see in our results. Far downstream 


(Fig. 6, right column, blue curves) the distribution is fea¬ 
tureless: the scattering is frequent enough that the distri¬ 
bution is isotropic. We emphasise that comparable pancake 
anisotropy in the downstream region is found for each Obu 
(Fig. 6) in the same momentum ranges in which a spike is 
seen; this suggests that such a local feature depends on local 
relative orientation of field line and shock surface and not 
by the value of the global obliquity Obu- 

Figure 6, left column, shows in detail the upstream 
anisotropy associated with the spike: particles spending a 
few gyroperiods in the thin shell ahead of the shock before 
being advected back to the shock or particles diffusing within 
the upstream shell and being accelerated. Upstream at the 
spike (—30r° < x < 0, Fig. 6, left column, red curves), we 
find two bumps at Uiocai — iO.45 for each Obu- In order to 
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interpret the bumps, we compare Fig. 7 with the adiabatic 
test-particle theory and previous simulations in weak tur¬ 
bulent field (Decker 1988). We should point out that this is 
expected to be a qualitative argument as the adiabatic test- 
particle theory is based on the conservation of the particle’s 
magnetic moment before and after a single encounter with a 
quasi-perpendicular shock whereas our simulations include 
both single and multiple encounters and track the pitch- 
angle also during the shock encounter. The lowest curve in 
Fig. 7 shows at a distance between [—450]r° and [—420]r° 
an excess at fiiocai = — 1 with an increasing anisotropy for 
smaller ; similar broadening was found in the comparison 
of the scatter-free limit with the case cr^ = 0.09 (see Fig. 17 
in Decker (1988))^. However, by sampling the distribution of 
IJiiocai closer to the shock (Fig. 7), we find a stronger depar¬ 
ture from the scatter-free limit. The bumps result from the 
turbulence that deforms the field lines locally at the shock 
front. Fluctuations allow a single field line to have multi¬ 
ple connection points along the shock front, leading to a 
bidirectional streaming and bumps in both hemispheres as 
reflected particles can come back to the shock. The bump at 
l^iocai < 0 is higher than fiiocai > 0 as a result of the excess of 
particles moving anti-parallel to the shock upstream in the 
scatter-free limit. In the upper curves (at [—30 : 0]r°) the 
bumps are enhanced and narrower for small cr^, although 
we point out that the fiiocai is greater than the limiting 
anisotropy in the scatter-free limit (/x < —0.81). Thus, Fig. 
6 shows that the presence of two peaks and f^iocai is inde¬ 
pendent on Obu- This corroborates our interpretation that 
the spike originates from the interaction of the small-scale 
magnetic field with the shock: in the presence of fluctuations 
the large-scale, or average, magnetic obliquity Obu is not as 
important close to the shock and the particle anisotropy is 
described by giiocai only. 

At higher momentum ([10 — 20]po), the thickness of the 
shell is too narrow to allow for pitch-angle isotropization: for 
instance, the red curves in Fig. 6 are calculated in a shell of 
thickness 30r°, i.e., only 1.5 times broader than the initial 
gyroradius for particles with momentum 20po, insufficient to 
provide particle scattering for isotropization even for a mean 
free path as small as the gyroradius, i.e., Bohm regime. 

The use of a pitch-angle with respect to the local field 
direction requires an additional justification. If we con¬ 
sider the parallel shock case (Fig. 4), the mean free path 
A|| ~ ^ 11/'^5 with V particle speed in the local plasma frame 
and ACII estimated from quasi-linear theory^, is larger or com¬ 
parable to spike thickness, so within the spike the particles 
cannot isotropize in the pitch-angle, calculated with respect 
to Bo, i.e. iJLgiohai — P • Bo/pHo. In other terms the steady 
state f{iigiobai) calculated close to the shock might reflect 
the isotropy in global used to inject particles at the shock 
rather than being a genuine isotropy produced by down¬ 
stream scattering: particles are not expected to behave dif¬ 
fusively near the shock, within the spike. In addition, we find 
that behind the shock f{figiobai) is flat (not shown in this 

For the parameters chosen here, at all energies, the adiabatic 
test-particle theory provides a strong anisotropy for upstream re¬ 
flected particles receding from the shock (/(p) non-vanishing for 
g < - 0 . 81 ) 

^ Same conclusion is found if ACy is numerically estimated from 
the exponential roll-over of the upstream particle density. 


Vq/Ui = 3.5, Gen = 85°, = 0.3 



P/Po 

Figure 8. Steady-state differential intensity spectrum as a func¬ 
tion of momentum in the local plasma frame in various spatial 
regions (within the two shells on both sides of the shock and 
far from it) for Ob^ = 85°, =0.3 and vq/Ui = 3.5 compared 

with DSA prediction far downstream. The top panel shows dJ/dE 
downstream (within 0 < x < 30r° the solid line and [0.1 : 
the dashed line). The bottom panel bear on the upstream flow 
(within —30r° < x < 0 the solid line and x < —30r° the dashed 
line). 

paper), for Obu = 0° at the same low momentum ([2 — 4]po) 
which exhibits the spike. However, in our simulations we 
compute the effect of a fluctuating magnetic field on a collec¬ 
tion of particles by summing up the instantaneous (on time- 
scale smaller than f^^^) position and pitch-angle of the ac¬ 
tual trajectories. This suggests that whether or not the dis¬ 
tribution isotropizes through scattering must be evaluated 
at a smaller scale with fiiocai = p-B/pH. Since we ensemble- 
average over several turbulence realizations and spatially av¬ 
erage over uncorrelated regions of the shock surface, any 
systematic effect of the turbulence on the anisotropy should 
average out. 


7 SPECTRA 

Figure 8 shows the steady-state differential intensity spec¬ 
tra dJ/dE = p^f{p) as a function of the momentum in the 
local plasma frame within various spatial regions compared 
with the downstream DSA prediction, i.e., test-particle limit 
power-law where a = 3r/(r—1) {a = 4.2 for r = 3.5). 

In the far downstream region the spectrum steepens to the 
DSA limit at smaller momentum for parallel than for per¬ 
pendicular case (see also Giacalone (2005)): since typically 
^ /^ll fewer energetic particles return to the shock, if 
highly oblique, to be further accelerated. A comparison of 
the spectra near the spike, but downstream of it, to that 
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far downstream of the shock is shown in the top panel of 
Fig. 8 within the range [3 — 20]po- There is a prominent 
MZocaZ-anisotropy in this case (see Fig. 6), leading to excess 
of low-energy particles at the spike with respect to the far 
downstream region®. In both cases considered, the spectrum 
at the spike softens because the population piling-up at the 
spike is deplenished in the most energetic particles. 


8 APPLICATION TO ASTROPHYSICAL 
SHOCKS 

In this Section, we discuss the implications of our findings 
for the in-situ measurements of interplanetary shocks^. 

8.1 Low-corona shock 

The detection of spikes at interplanetary shocks has been 
rare so far, due to limited time-resolution of the energetic 
particles detectors. We predict in Fig. 9 the energetic parti¬ 
cles profiles at travelling shocks moving at high speed close 
to the Sun. For protons injected at kinetic energy = 50 
keV in the slow solar wind at 10 solar radii {Bq = 0.025 
G) at a relatively strong (r = 2.5) quasi-perpendicular 
— 85°) shock we use Ui = 2,500 km s“^ (Ma = 2) 
and Lc = 10“^ AU. The predicted duration of the shock 
spike as seen by a spacecraft crossing the shock along the 
x-axis turns out to be a few seconds (shorter for greater Obu 
as shown previously). Finally, we find comparable amplitude 
and duration of the spike if the momentum is calculated in 
the local plasma frame or in the shock frame. Instruments 
on the upcoming NASA mission Solar Probe Plus, includ¬ 
ing those that measure energetic particles (McComas et al. 
2014), will have a cadence less than this and will be capable 
of resolving the shock spike events. 

8.2 Solar Termination Shock 

The crossing of the solar termination shock (TS) by 
Voyager-1 (Decker et al. 2005) reveals that energetic par¬ 
ticles exhibited a spike in intensity at the time of shock 
passage in high-energy protons (3.4 — 17.6 MeV) and ions 
(40 — 53 keV). In Fig. 10 we display the simulated density 
profile with parameters compatible to the numerical solu¬ 
tion of the FTE for Obu = 80° in Florinski et al. (2008): 
Ui = 350 km s“^, r = 3.0, vo/UicosObu = v^5/3; also, here 

® Such an excess can be quantified with a non-linear least-squares 
fit: for the case 9Bn = 85°, = 0.3 (Fig. 8) we find slopes 

— 1.393T0.045 (spike downstream) and —1.294A0.044 (far down¬ 
stream), compared with the DSA value —2.2 expected at higher 
momenta. For the case Obu = 0°, cr^ = 1.0, we find slopes 

— 1.680±0.022 (spike downstream) and —1.573 ±0.055 (far down¬ 
stream); clearly, the parallel shock is closer to DSA limit. 

^ We have performed simulations of supra-thermal protons at 
shocks of galactic supernova remnant as well. We find a shock 
spike for injection energy E = 10 MeV at a strong shock (r = 4) 
propagating into the ISM (Bq = SgG with small fluctuation, i.e., 
(7^ = 0.3, Lc = 2.2 X 10^^ cm) with Obu = 85° and Ui = 6,600 
km/s. However, the thickness of the spike is very narrow (10“® 
pc) compared with the typical radius of a few parsec for a young 
or middle-aged remnant. 


Vq/Ui = 3.9, 0Bn = 85, =1.0, r = 2.5 

Time (s) in spacecraft frame 



-600 -400 -200 0 200 400 600 


Figure 9. Steady-state spatial density profiles of accelerated par¬ 
ticles for various particle momentum ranges at = 85° and 
(7^ = 1.0. The density is normalised to the far downstream limit. 
In the lower x-axis the distance is in units of initial upstream gy- 
roradius, in the upper the time is in units of a spacecraft past by 
the shock. The spike duration between 1.5 and 8 MeV amounts to 
(~ 1 sec) ensuring observability by an energetic particle detector 
with current and future time-resolution. 


Bo = 5 X 10“^ G, Ma — 10 and Lc — 0.01 AU. Ions in 
108 — 180 keV (lower panel in Fig. 10) exhibit a shock spike 
less prominent than the spike at 150 keV found in Florinski 
et al. (2008), i.e., a factor 10 (a smaller spike is found 
for smaller Obu)- The finite energy interval chosen here, 
in contrast with a mono-energetic population in Florinski 
et al. (2008), might increase the far downstream limit rel¬ 
atively to the spike, reducing the amplitude of the spike. 
However, we have used a 3D-isotropic magnetic turbulence 
power spectrum, in contrast with the slab turbulence used 
in Florinski et al. (2008): as a general consequence of par¬ 
ticle transport (Jokipii et al. 1993), in a slab turbulence 
particles are unphysically forced to adhere to the field lines, 
i.e., to a direction lying on the shock surface (in the case of 
Obu — 80° considered in Florinski et al. (2008)), on a time- 
scale shorter than the scattering time (Fraschetti Jokipii 
2011). Thus, one would expect a less prominent spike in 
the slab case as the advection carries particles away from 
the shock. Although a slab might correctly approximate the 
Alfvenic turbulence at the TS, the isotropic turbulence used 
here allows to account for both perpendicular and cross field 
diffusion with no assumptions. 
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Vq/Ui =7.6, 0Bn = 8O°, o^ = 0.1 
AU 



-800 -600 -400 -200 0 200 400 600 800 


Figure 10. Spatial density profiles of solar termination shock 
accelerated particles for various particle kinetic energy ranges at 

— ^0° and (7^ = 0.1. The density is normalised to the far 
downstream limit. The distance is in the lower x-axis in units of 
initial upstream gyroradius, in the upper in units of AU. 

9 DISCUSSION AND CONCLUSION 

We have performed several test-particle numerical simu¬ 
lations by integrating the trajectories of individual ener¬ 
getic particles scattering off a magnetic turbulence advected 
through a planar shock to study the formation of spike¬ 
like structures, or localised intensity enhancements, not pre¬ 
dicted by the standard steady-state DSA. Previous works 
based on numerical integration of the focused transport 
equation (le Roux et al. 2007) argued that no local enhance¬ 
ments can be found at parallel, or small obliquity, shocks for 
low-energy suprathermal particles. Our direct particle tra¬ 
jectory integration shows that such spikes do indeed exist 
at quasi-parallel shocks that move through large-scale mag¬ 
netic turbulence. 

The presumed continuous transition from SDA to DSA 
as dominant scenario of particle acceleration at shocks as 
going from quasi-perpendicular to quasi-parallel geometry, 
although roughly consistent with several observations and 
models, contains discrepancies and remains to be under¬ 
stood (see Lee et al. (2012) for a review in the context 
of interplanetary shocks). Our results aim at reconciling 
the dichotomy between the modelling of spike at quasi¬ 
perpendicular shocks with SDA (Decker 1983), and the ar¬ 
gument for less oblique shocks {0 < 70°) that the shock en¬ 
hancement arises from particles trapped in the upstream by 
self-generated hydromagnetic waves, with a resulting drop 
of particles mean free path (Scholer 1985). 

We find a local pitch-angle anisotropy resulting in a 
pancake-like shape, independent on the value of the large- 
scale obliquity Obu- Therefore it seems not necessary to as¬ 
sume the presence of self-generated waves as the turbulence, 
inherent property of the medium travelled by shocks, is able 


to trap particles also at Obu ~ 0° shocks, at least at supra¬ 
thermal energies. 

The approach of simulating the direct particle trajec¬ 
tory in magnetic turbulence used here naturally includes 
parallel and perpendicular diffusion to the large scale uni¬ 
form magnetic field, without assuming a particular form of 
perpendicular transport. We note that a decorrelation of 
low-energy particles might occur within a few gyroperiods 
time-scale even in weak turbulence {{5B/BqY — 0.1), as di¬ 
rect simulations in Fraschetti & Giacalone (2012) suggest. 
Suppressing the perpendicular diffusion, namely artificially 
constraining the particle motion (see for instance Fraschetti 
& Jokipii (2011)), might unphysically limit the spike ampli¬ 
tude in the quasi-perpendicular case because the advection 
along the flow artificially dominates over the perpendicular 
diffusion (Giacalone & Jokipii 1999). 

The spatial features at proton-scale sought here do not 
need hybrid simulations as test-particle approach can con¬ 
sistently capture the phenomenon we are interested in. How¬ 
ever, in test-particle simulations the definition of magnetic 
obliquity has only large-scale validity, because even at pro¬ 
ton gyroscales the real shock structure will be corrugated 
by plasma micro-instabilities, so that the same particle re¬ 
flected to the upstream will experience a different Obu at its 
successive encounters with the shock. The shock planarity 
is also affected by the density of the medium it is travel¬ 
ling into: large-scale corrugation induced by the upstream 
inhomogenieties will also convert bulk kinetic energy of the 
shock to magnetic fluctuation power at the energetic proton 
gyroscale, reducing energetic particles mean free path, as 
demonstrated for high-Mach number non-relativistic shocks 
(Giacalone & Jokipii 2007; Fraschetti 2013) and for mildly 
relativistic shocks (Mizuno et al. 2014). Gross-shock electro¬ 
static potential, generated by the charge separation between 
ions and electrons at the shock due to differential deflection 
in the magnetic field, has been neglected in this paper as 
confined to the electron skin-depth scale. Thus, the consis¬ 
tent study of electron density spatial profile, deferred to a 
future work, requires the use of hybrid simulations (particle 
electrons and fluid ions). 

We have neglected the effect on the spike of the up¬ 
stream self-generated waves. The amplitude of the self¬ 
generated turbulence grows from perpendicular to parallel 
shock (Lee et al. 2012). Therefore, should spike-events down¬ 
stream in parallel shocks originate upstream from hydro- 
magnetic waves trapping and accelerating particles which 
are further convected downstream and produce spikes, one 
should find the largest amplitude spikes at parallel shocks. 
However, in our simulations spikes at globally oblique shocks 
{Obu — 70°) reach the highest amplitude. This feature, not 
explained in our test-particle approximation, requires fur¬ 
ther investigation. 

In summary, since spikes are shown here to be an in¬ 
herent property of shocks, we underline that, beside the 
Voyager-1 crossing of termination shock, a much larger sam¬ 
ple is expected to be seen at interplanetary shocks by high¬ 
time resolution particle detectors. To this purpose, investi¬ 
gation of time-dependent spatial profiles, not considered in 
this work, might also help interpretation of data. Our simu¬ 
lations might be relevant to the future Solar Probe Plus and 
Solar Orbiter Missions. Unprecedently high time-resolution 
in situ measurements in the solar atmosphere (high solar 
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corona down to ~ 10 solar radii) will enable measurements 
of small-scale spatial features in energetic electrons and pro¬ 
tons distributions at interplanetary shocks. This will con¬ 
tribute to clarify the behaviour of non-thermal particles 
at non-relativistic shocks, one of the main focuses of So¬ 
lar Probe Mission, possibly impacting on the understanding 
of energetic particles at supernova remnant shocks. 


ACKNOWLEDGMENTS 

We acknowledge the useful discussions with J. R. Jokipii and 
J. Kota. We thank the anonymous referees for useful sug¬ 
gestions and comments. This work was supported, in part, 
by NASA under Grants NNX10AF24G and NNX11A064G. 
Work by JG was also supported, in part, from the ISIS in¬ 
strument on NASA’s Solar Probe Plus mission, and by the 
NSF under grants AGS 1135432 and AGS 1154223. This work 
benehted from technical support by the computer cluster 
team at the Department of Planetary Sciences at University 
of Arizona. 


REFERENCES 

Acciari V. A. et ah, 2011, ApJ, 730, L20 
Albert J., et ah, 2007, ApJ, 664, L87 
Armstrong T. P., Krimigis S. M., Behannon K. W., 1970, 
J. Geophys. Res., 75, 5980 

Axford W. I., Reid G. G., 1963, J. Geophys. Res., 68, 1793 
Bell A. R., 1978, MNRAS, 182, 147 
Decker R. B., 1983, J. Geophys. Res., 88, 9959 
Decker R. B., 1988, Space Science Reviews, 48, 195 
Decker R. B., et al. 2005, Science, 309, 2020 
Drury L. O’ G., 1983, Rep. Prog. Phys., 46, 973 
Earl J. A., Jokipii J. R., Morhll J., 1988, ApJ, 331, L91 
Ellison D. G., Baring M. G., Jones E. G., 1995, ApJ, 453, 
873 

Ellison D. G., Baring M. G., Jones, E. G., 1996, ApJ, 473, 
1029 

Erdos G., Balogh A., 1994, ApJS, 90, 553 
Elorinski V., Zank G., le Roux J. A., 2008, Adv. Space 
Res., 41, 361 

Eraschetti E., 2013, ApJ, 770, 84 
Eraschetti E., Giacalone, J., 2012, ApJ, 755, 114 
Eraschetti E., Jokipii J. R., 2011, ApJ, 734, 83 
Giacalone J., 2005, ApJ, 624, 765 
Giacalone J., 2012, ApJ, 761, 28 
Giacalone J., Jokipii J. R., 1999, ApJ, 520, 204 
Giacalone J., Jokipii J. R., 2007, ApJ, 663, L41 
Giacalone J., Jokipii J. R., 2009, ApJ, 701, 1865 
Gieseler U.D.J., Kirk J.G., Gallant Y.A., Achterberg, A., 
1999, A&A, 345, 298 

Isenberg P. A., 1997, J. Geophys. Res., 102, 4719 
Jokipii J. R., 1982, ApJ, 255, 716 
Jokipii J. R., Giacalone J., 2007, ApJ, 660 ,336 
Jokipii J. R., Kota J., Giacalone J., 1993, Geophys. Res. 
Lett., 20, 1759 

Jones E. G., Birmingham T. J., Kaiser T. B., 1978, Phys. 
Eluids, 21, 347 


Kota J., Jokipii J. R., 2004, in AIP Gonf. Proc. 719, Physics 
of the Outer Heliosphere, ed.V. Elorinski,N.V. Pogorelov, 
& G. P. Zank (Melville:AIP), 272 
Lario D., Ho G. G., Decker R. B., Roelof E. G., Desai M. L, 
Smith G. W., 2003, in AIP Gonf. Proc. 679, Solar Wind 
Ten, ed. M. Velli, R. Bruno, & E. Malara, 640 
le Roux J. A., Webb G. M., Elorinski V., Zank G. P., 2007, 
ApJ, 662, 350 

Lee, M. A., Mewaldt, R. A., Giacalone, J. 2012 Space Sci. 
Rev., 173, 247 

McGomas et ah, 2014, Space Science Reviews, in press 
Medvedev M. V., Loeb A., 1999, ApJ, 526, 697 
Mizuno Y., Pohl M., Niemiec J., Zhang B., Nishikawa K.-L, 
Hardee P. E. 2014, MNRAS, 439, 3490 
Naito T., Takahara E., 1995, MNRAS, 275, 1077 
Ostrowski M., 1991, MNRAS, 249, 551 
Parker E. N., 1958, ApJ, 128, 664 
Parker E. N., 1965, P&SS, 13, 9 

Press W. H., Elannery B. P., Teukolsky S. A., Vetterling 
W. T., 1986, Numerical Recipes (Gambridge: Gambridge 
Univ. Press) 

Schlickeiser, R. & Vainio, R., 1999, Astrophysics and Space 
Science, 264, 457 

Scholer M., 1985, Gollisionless Shocks in the Heliosphere: 

Reviews of Gurrent Research, 287 
Singer S., Montgomery M. D., 1971, J. Geophys. Res., 76, 
6628 

Zuo P.-B., Peng X.-S., 2013 Ghinese Phys. Lett. 30 019601 


APPENDIX A: PARTICLE TRAJECTORY 
INTEGRATION 


The equation of motion of a test-particle with charge e and 
mass m moving with velocity v{t) = | v(t) | in a magnetic held 
B(x) is the Lorentz equation, written in the shock frame as 


I u(i) X n(x,i) 

dt 


(Al) 


where we dehned D(x, t) = eB(x, t)/(?Tic) with the Lorentz 
factor 7 (t) = 1/y^l — (v{t)/c)^ , c the speed of light in vac¬ 
uum and u{t) = 7 v(t)/c. The local electric held E(x, t) on 
the z-th side of the shock is calculated within the ideal in- 
hnitely conductive MHD approximation in the shock frame: 
^(x, t) = eE(x, t)/(mc) = — (Ui)/c x D(x, t). We calculate 
the trajectory of the particle in a magnetic held as a so¬ 
lution of Eq.(Al) where t is the time in the rest frame of 
the plasma, coincident with time in the shock frame be¬ 
cause Ui c. The quantity D(x, t) in Eq.(Al) is given by 
D(x, t) = Do + JD(x, t) where Do = (e/mc7)Bo, in terms of 
the background magnetic held Bo constant and statistically 
uniform, and (5D(x, t) is the turbulent component varying in 
the three-dimensional space; the time variation is accounted 
for by the frozen-in condition. 

We use the prohle of the bulk velocity specihed in 
Sect. 3, namely Ux{x) = Ui for x < 0, Ux{x) = U 2 for x > 0. 
The magnetic held embedded in the plasma is solution of 
the induction equation in the ideal MHD approximation, ad- 
vected along the how (Jokipii & Giacalone 2007; Giacalone 
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& Jokipii 2009): 


Bx(x,y,z,t) = Bx(xo,y,z,to) 

By(x,y,z,t) = By{xo,y,z,to)Ui/U(x) 

BFx,y,z,t) = BFxo,y,z,to)Ui/U{x) (A2) 

where xq and to are related by the characteristics equation: 
t — to = dx' jU{x'). For convenience, we choose xq in the 

far upstream of the shock and use the characteristic equation 
to find to. If B is known at {xo,y, z,to), using A2 the field 
can be calculated at every point where B is advected: at any 
time the magnetic field at particle position depends on the 
previous history of the field line carried by the flow. 

The fluctuating field comprises of transverse waves at 
every point of physical space occupied by the particle, each 
with a random amplitude, phase, orientation and polar¬ 
ization (Giacalone & Jokipii 1999; Fraschetti & Giacalone 
2012). We use a three-dimensional isotropic, or Kolmogorov, 
power spectrum in the upstream plasma. The coherence 
scale Lc of the turbulence is chosen such that rgjLc ^ 1. 
The discrete modes are logarithmically equispaced. As in 
Fraschetti & Giacalone (2012) the randomization is per¬ 
formed over the initial particle velocity orientation in the 
local plasma frame and over the turbulence realization, shuf¬ 
fled every 50 particles. Including the fluctuating field statis¬ 
tics ensures a meaningful comparison with a theoretically 
computed turbulence power spectrum, which is by defini¬ 
tion an average over an ensemble of field realizations. 

We inject a large number of particles on the plane 
X = 0.5r°, with randomly oriented velocity, but fixed plasma 
frame initial kinetic energy. We integrate the Eq. (Al) by 
using a time-step adjustable Burlisch-Stoer method (Press 
et al. 1986). The integration of the equation of motion is per¬ 
formed in the shock frame which allows to calculate explic¬ 
itly at every time step the ideal MHD electric field specified 
above. This approach is well suited to investigate the particle 
transport in turbulence. A different approach of computing 
different realizations of the turbulent magnetic field in ev¬ 
ery cell of an Eulerian grid is much more time-consuming. 
The latter approach would also require adapting the lattice 
spacing in order to maintain the same space resolution in 
physical space. Also, we did not implement the “particle¬ 
splitting” method as the momentum boundary is relatively 
small {p < Pb = 500po). 

In the steady state the boundary conditions are to be 
given in position and momentum (Giacalone 2005): the tra¬ 
jectory of each particle is calculated until it reaches a far 
downstream return spatial boundary at Xb = 2.5x lO^f/i/Qo, 
or until the momentum p in the plasma frame becomes 
greater than the upper boundary pb = 500po. If P > Pb, the 
particle is removed from the system. If x > we use the 
probability return algorithm (Ellison et al. 1996; Giacalone 
2005): if the particle passes the test, it is reinitialised at 
the boundary plane x = Xb with the same kinetic energy 
and same p, z coordinates held at the last time-step before 
the boundary crossing with an isotropic pitch-angle distribu¬ 
tion in the local downstream plasma frame. The re-injection 
mimics the returning flux of energetic particles at the loca¬ 
tion Xb from the far downstream (Ellison et al. 1996). In fact, 
if we replace the downstream boundary at x = x?, with an 
open boundary, at a location close to the boundary x < Xb 
we find the expected isotropic distribution. In these numeri¬ 
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Figure Al. Numerical convergence for various sampling interval 
Ax/Vg = 2, 4,12 in all momentum bands considered in this paper. 

The small panels on the left zoom into the shock in the same units. 


cal experiments we used Xb > Lc- This is consistent with the 
fact that the coordinates (po,^o) in the injection plane are 
spread over a scale greater than Lc'. the simulation samples 
various regions of the shock both in the shock plane and in 
the perpendicular direction. 

The steady-state density profiles comprise pseudo¬ 
particles: with a fixed cadence in time (AtQo = 5) and given 
a particle momentum range in the local plasma frame, we 
sample the position of each particle in spatial bins with uni¬ 
form thickness along the x-axis from —Xb to Xb until each 
particle is removed from the system (the resulting spikes in 
the density profiles are verified to be roughly independent 
on the thickness of the bins). The density profile within the 
given momentum range is found by summing up in every 
spatial bin the total number of pseudo-particles within the 
given momentum range collected at every time. The inter¬ 
val AtLlo is conveniently chosen not too large to track the 
gyro motion of the particles around the shock and not too 
small to reduce computational time. Likewise, the spectrum 
is produced by sampling the momentum of each particle in 
momentum bins with uniform thickness in logarithmic scale 
within [0.9po : Pb] with cadence AtLlo = 5. Eigure Al shows 
the numerical convergence of the spatial profiles for various 
spatial sampling interval Ax/r°. 
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